/* libwatopt++ - Optimization Toolkit
 *
 * Copyright (c) 2007 Christopher Alexander Watford
 * <christopher.watford@gmail.com>
 *
 * Permission is hereby granted, free of charge, to any person obtaining a copy
 * of this software and associated documentation files (the "Software"), to
 * deal in the Software without restriction, including without limitation the
 * rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
 * sell copies of the Software, and to permit persons to whom the Software is
 * furnished to do so, subject to the following conditions:
 *
 * The above copyright notice and this permission notice shall be included in
 * all copies or substantial portions of the Software.
 *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
 * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
 * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
 * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
 * IN THE SOFTWARE.
 *
 * $Id: ga.cpp 5 2007-05-12 19:19:02Z christopher.watford $
 */

#include <cstdlib>
#include <algorithm>
#include <functional>
#include <numeric>
#include <vector>
#include <cmath>

#include <opt.hpp>
#include <ga.hpp>
#include <utils.hpp>

namespace watopt
{

   GeneticAlgorithmOpt::GeneticAlgorithmOpt(size_t populationSize,
      size_t generations)
   : _m_PopSize(populationSize), _m_MaxGenerations(generations)
   {
   }

   void
   GeneticAlgorithmOpt::Reaper(std::vector<IFitnessImplementation*>& pool)
   {
      while(pool.size() > _m_PopSize)
      {
         IFitnessImplementation* member = pool.back();
         pool.pop_back();
         
         delete member;
      }
   }

   bool
   GeneticAlgorithmOpt::EndBreeding(size_t breedings, size_t children)
   {
      return (children < _m_PopSize/2 ? false : true);
   }

   void
   GeneticAlgorithmOpt::Selection(
      const std::vector<IFitnessImplementation*>& pool,
	   IFitnessImplementation** parent1, IFitnessImplementation** parent2,
      double minFit, double maxFit, double avgFit, double variance)
   {
	   // select 2 parents by a tournament selection
	   if(parent1 && parent2)
	   {
         std::vector<IFitnessImplementation*> parents;

         // 1. pick 10% of the population at random to be parents
         size_t poolSize = pool.size();
         size_t parentCount = std::max<size_t>(poolSize/10,2);
         while(parents.size() < parentCount)
         {
            size_t p1 = rand_index(poolSize);
            while( std::find(parents.begin(), parents.end(),
                              pool[p1]) == parents.end() )
               p1 = rand_index(poolSize);

            parents.push_back(pool[p1]);
         }

         // 2. sort them by fitness
		   std::sort(parents.begin(), parents.end(),
                     IFitnessImplementation::FitnessPtrComp());

         // 3. take the two best parents
		   (*parent1) = parents[0];
		   (*parent2) = parents[1];
	   }
   }

   static double
   fitAccumulate(double fit, IFitnessImplementation *ff)
   {
	   return fit + (ff ? ff->Fitness() : 0.0);
   }

   /**
    * Default populace convergence function does nothing. Suggested
    * implementations would be to have a natural disaster occur which would
    * reap the population down to a few of the high performing members
    * and then reinsert random members (perhaps running a local search on
    * them).
    *
    * \param pool Population which has converged
    */
   bool
   GeneticAlgorithmOpt::PopulaceConverged(
      std::vector<IFitnessImplementation*>& pool)
   {
      return false;
   }

   /**
    * Initialize the algorithm by filling the population with random
    * members. Then sort them by fitness.
    */
   void
   GeneticAlgorithmOpt::Initialize(void)
   {
	   while(_m_Population.size() < _m_PopSize)
		   _m_Population.push_back(Generate());

      std::sort(_m_Population.begin(), _m_Population.end(),
                  IFitnessImplementation::FitnessPtrComp());

      GatherStatistics(_m_Population);
   }

   /**
    * Compute statistics on the population. Does minimum fitness, maximum
    * fitness, average fitness, and variance.
    *
    * \param pool Population to gather statistics on
    */
   void
   GeneticAlgorithmOpt::GatherStatistics(
      const std::vector<IFitnessImplementation*>& pool )
   {
      _m_MinFit = _m_Population[0]->Fitness();
      _m_MaxFit = _m_Population[_m_Population.size()-1]->Fitness();
	   _m_AvgFit = std::accumulate(_m_Population.begin(), _m_Population.end(),
                     0.0, fitAccumulate) / _m_Population.size();

      _m_Variance = 0.0;
      std::vector<IFitnessImplementation*>::const_iterator it;
      for(it = pool.begin(); it != pool.end(); ++it)
      {
         double fitdelt = (*it)->Fitness() - _m_AvgFit;
         _m_Variance +=  fitdelt*fitdelt;
      }
   }

   /**
    * Basic GA solver.
    *
    * 1). Select parents via tournament selection
    * 2). Crossover children and add them to a secondary population
    * 3). Repeate 1 and 2  until the population size increases by 50%
    * 4). Return children to the population
    * 5). Mutate the population
    * 6). Sort the population by fitness
    * 7). Trim the population down to its prebreeding size
    * 8). Gather statistics on the population
    * 9). Repeat 1 through 8 until Max Generations is reached
    * 10). Return the top population member as the solution
    */
   IFitnessImplementation*
   GeneticAlgorithmOpt::Solve(void)
   {
      if(_m_Instrument)
         _m_Instrument->Start(_m_Population, _m_MinFit, _m_MaxFit, _m_AvgFit,
                                 _m_Variance);

	   for(size_t GG = 0; GG < _m_MaxGenerations; GG++)
	   {
		   std::vector<IFitnessImplementation*> newPop;

		   // 1. Breed population
		   for(size_t bb = 0; EndBreeding(bb, newPop.size()); bb++)
		   {
			   IFitnessImplementation *p1, *p2;
   			
			   // 2a. Select 2 parents
			   Selection(_m_Population, &p1, &p2, _m_MinFit, _m_MaxFit, _m_AvgFit,
                        _m_Variance);

            if(_m_Instrument)
               _m_Instrument->ParentsSelected(p1, p2);

			   // 2b. Create children, add them to the new population
			   Crossover(p1, p2, newPop);
		   }

		   // 2. Reintroduce children into population
		   _m_Population.insert(_m_Population.end(), newPop.begin(),
                                 newPop.end());
		   newPop.clear();

		   // 3a. Mutate
		   std::for_each(_m_Population.begin(), _m_Population.end(),
            std::bind1st(
               std::const_mem_fun1_t<void,
                                     const GeneticAlgorithmOpt,
                                     IFitnessImplementation*>(
                                       &GeneticAlgorithmOpt::Mutate
               ),
               this
            )
         );

		   // 3b. Resort
		   std::sort(_m_Population.begin(), _m_Population.end(),
                     IFitnessImplementation::FitnessPtrComp());

		   // 4. Reap population down to size
		   Reaper(_m_Population);

		   // 5. Find min, average, and max fitness
         GatherStatistics(_m_Population);

         // 6. Determine if the population has converged
         if(_m_MinFit == _m_MaxFit)
         {
            if(_m_Instrument)
               _m_Instrument->PopulaceConverged(GG, _m_Population, _m_MinFit);

            // 6a. Run the population converged function and if it returns
            // true, gather statistics again because it altered the
            // population
            if(PopulaceConverged(_m_Population))
               GatherStatistics(_m_Population);
         }

         if(_m_Instrument)
            _m_Instrument->EndGeneration(GG, _m_Population, _m_MinFit,
                                          _m_MaxFit, _m_AvgFit, _m_Variance);
	   }

      if(_m_Instrument)
         _m_Instrument->Stop(_m_Population, _m_MinFit, _m_MaxFit, _m_AvgFit,
                              _m_Variance);

	   return _m_Population[0];
   }

   /**
    * Get all of the population members.
    *
    * \param solutions Vector to fill with all of the solution members
    * does not clear the input vector.
    */
   void
   GeneticAlgorithmOpt::GetSolutions(
      std::vector<IFitnessImplementation*>& solutions ) const
   {
      solutions.insert(solutions.end(), _m_Population.begin(),
                        _m_Population.end());
   }
}
